

GetK = function(a, l, ratio, k.max, K){

    # one ratio for multiple a, l
        
    r  = ratio["l0",1]^l[,"0"] * ratio["ak",l[,"0"]+1]^a[,"0"] 
    for(k in 2:K){
        r = r * ratio["lk", a[,k-1] + 1]^l[,k] *
            ratio["ak", l[,k] + 1]^a[,k]
    }
    r = r * ratio["lK", a[,K] + 1]^l[,K+1] *
        ratio["aK", l[,K+1] + 1]^a[,K+1]
    
    k = r/k.max
    
    return(k)
}

GetKGen = function(a, l, ratio, k.max, K){

    # multiple ratio, each for one a, l

    r  = ratio[,"l0",1]^l[,"0"] * diag(ratio[,"ak",l[,"0"]+1])^a[,"0"]
    for(k in 2:K){
        r = r * diag(ratio[,"lk", a[,k-1] + 1])^l[,k] *
            diag(ratio[,"ak", l[,k] + 1])^a[,k]
    }
    r = r * diag(ratio[,"lK", a[,K] + 1])^l[,K+1] *
        diag(ratio[,"aK", l[,K+1] + 1])^a[,K+1]

    k = r/k.max

    return(k)
}


power = function(x,y) x^y

GetKGen2 = function(a, l, ratio, k.max, K){
    
    # multiple ratio, for the same multiple a, l
    
    r  = t(outer(ratio[,"l0",1],l[,"0"],power)) * 
        t(ratio[,"ak",l[,"0"]+1])^a[,"0"] 
    for(k in 2:K){
        r = r * t(ratio[,"lk", a[,k-1] + 1])^l[,k] *
            t(ratio[,"ak", l[,k] + 1])^a[,k]
    }
    r = r * t(ratio[,"lK", a[,K] + 1])^l[,K+1] *
        t(ratio[,"aK", l[,K+1] + 1])^a[,K+1]
    
    k = t(r)/k.max
    
    return(k)
}


